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ABSTRACT 

Wc estimate the main-sequence age and heavy-element abundance of the Sun by 
means of an asteroseismic calibration of theoretical solar models using only low- 
degree acoustic modes from the BiSON. The method can therefore be applied also 
to other solar-type stars, such as those observed by the NASA satellite Kepler and 
the planned ground-based Danish-led SONG network. The age, 4.60±0.04 Gy, ob- 
tained with this new seismic method, is similar to, although somewhat greater than, 
today's commonly adopted values, and the surface heavy-element abundance by mass, 
Zs=0.0142±0.0005, lies between the values quoted recently by Asplund et al. (2009) 
and by Caffau et al. (2009). We stress that our best-fitting model is not a seismic model, 
but a theoretically evolved model of the Sun constructed with 'standard' physics and 
calibrated against helioseismic data. 
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1 INTRODUCTION 

The only way by which the age of the Sun can be estimated 
directly to a useful degree of precision is by accepting the 
basic tenets of solar-evolution theory and measuring those 
aspects of the structure of the Sun that are predicted by 
the theory to be indicators of age. We recognize that there 
are also indirect methods based on the more reliable dcter- 
mination of the age s of meteorites (e.g. |A mclin et al. 200^; 

■ ■ ■ M 



I Jacobsen et al ] |2008l : iBouvier fc WadhwaTBoiOf ). We recog- 
nize also that there is not a precise origin of time, such as 
a moment that one can uniquely define to be the time at 
which the Sun arrived on the main sequence. However, after 
initial transients, the central hydrog en abundance Xc de- 
clined almost linearly with time (e.g. lGoughlll995l ). so one 
can extrapolate Xc{t) backwards quite well to the time when 
Xc = Xo, the initial hydrogen abundance. That is the time 
that we adopt as our fiducial origin. A potential goal of fu- 
ture investigations of the type we describe here could be 
to ascertain whether the Sun arrived on the main sequence 
before the rest of the solar system formed, or at the same 
time. Unfortunately we have not yet succeeded in resolv- 
ing the matter, partly because the data errors are not yet 
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small enough, but mainly, as we discuss in § 4, because the 
uncertainties in the modelling are too great. 

The solar structure measurements must be carried out 
seismologically, and one is likely to expect greatest reliability 
of the results when all the available pertinent helioseismic 
data are employed. Of these, the most pertinent are the fre- 
quencies of the modes of lowest degree, because it is they 
that penetrate the most deeply into the energy-generating 
core where the helium-abundance variation records the in- 
tegrated history of nuclear transmutation. Moreover, it is 
also only they that can be measured in other stars. There- 
fore, there has been some interest in calibrating theoretical 
stellar models using only low-degree modes - here we use 
modes of degrees l—Q , 1, 2 and 3. The pros p ect was firs t 
discussed in detail bv [Christensen-DalsgaardI () 19841 . [TqSSI I. 
lUlrichI (| 19861 ) and ICoughl (|l987f r although prior to that it 
had already been pointed out that the helioseismic frequency 
data that were available at the time indicated that either 
the initial helium abundance Yq, or the age tQ, or both, 
are somewha t greater than the generally accepte d values 
(|Gough|[l983 : see also iGough fc Kosovichevlll990l). Subse- 
quent , more careful, calibrat i ons w er e discussed bvlCruenther 
(1989 [ ) iGough fc Novotn^ (Il990tl. iGuenther fc Demargue 
( 19971) IWeiss fc Schlattll (ll998D,lDziembowski et all (ll999D . 
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the la st have been reviewed by IChristensen-DalsgaardI 
l|2009l ). who dicusses some of the obstacles that need to be 
surmounted. Most of the calibrations did not address the 
influence of uncertainties in chem ical composition on the 
determination of tQ ; for example, IWeiss fc Schlattl (| 19981 ) 
adopted in their calibration the helioseismically determined 
values for the helium abundance in the convection zone, to- 
gether with the convection- zone depth. 

As a main-sequence star ages, helium is produced in 
the core, increasing the mean molecular mass fj, preferen- 
tially near the centre, and thereby inducing a local positive 
gradient of the sound speed. The resulting functional form 
of the sound speed c(r) depends not only on age Iq but also 
on the relative augmentation of fj,{r), which itself depends 
on the initial absolute value of n, and hence on the chemical 
composition: directly on the initial helium abundance Yq, 
via the equation of state, and, to a lesser degree, Zq, and 
indirectly, via the model calibration to the observed values 
R© and Lq of the radius R and the lum inosity L, on Zq 
and, to a lesser degree, Yn. iGoughl (|200ll ) tried to separate 
these two dependencies using the degree dependence of the 
small separation d„j = 3{2l + (vnj — i'„-i,i+2) between 
cyclic multiplet frequencies i'n,i, where n is order and I is de- 
gree. This is possible, in principle, because modes of differ- 
ent degree and similar frequency sample the core differently. 
However, the difference between the effects of Iq and Yq on 
the functional form of c(r) in the core is not very great, and 
consequently the error in the calibration produced by errors 
in the observed frequency data is uncomfortably high, as 
is also the case when a mean value of the large separation 
{i^n.i — Vn- i.i) is used in conjunctio n with the mean small 
separation l|Gough fc Novotnvlll990l ). 

This lack of sensitivity can be overcome by using, in 
addition to core-sensitive seismic signatures, the relatively 
small oscillatory component of the eigenfrequencies induced 
by the sound- speed glitch associated with helium ionization 
l|Goughl bOOj ) , whose amplitude is close to being propor- 
tional to helium abundance Y jHoudek fc Goughl l2007lj ) . 
The neglect of that component in the previously employed 
asymptotic signature had not only omitted an important 
diagnostic of Y , but had appeared to imprint an oscilla- 
tory contamination in the calibration as the limits (fci,A;2), 
whe re fc = n -\- ^l, of the adopted mode range was var- 
ied llGoughll200ll ). It therefore behoves us to decontaminate 
the core signature from glitch contributions produced in the 
outer layers of the star (from both helium ionization and the 
abrupt variation at the base of the convection zone, and also 
from hydrogen ionization and the superadiabatic convective 
boundary layer immediately beneath the photosphere). To 
thi s end a helioseism i c glitch signature has been developed 
bv lHoudek fc Goughl (|2007bt ). from which its contributions 
5vn,i to the frequencies can be computed and subtracted 
from the raw frequencies Vn,i to produce effective glitch-free 
frequencies Vsn,i to which a glitch-free asymptotic formula - 
equation - can be fitte d. The solar calibration is then ac- 
complished as previously (|Goughll200ll ) by fitting theoretical 
seismic signatures to the observations by Newton-Raphson 
iteration, using a carefully computed grid of calibrated mod- 
els to compute derivatives with respect to Zq and the age 
of each model. The result of the first preliminary calibration 
by this method, using Bi SON data, has been reported by 
iHoudek fc Goughl (|2007al ). Here we enlarge on our discus- 
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Figure 1. Functional forms fx of the integrands ipx in X = 
(t>x'ir, where fxii") = 'f>x (r)/\4'x(0)\ and where X = A, C , 
or F, plotted for Model S of IChristensen-Dalsgaard et al. I lll996l ') 
over the inner half of the interval (0, R) of r. The parameters A, 
C and F are sensitive particularly to the structure of the core, 
being progressively more centrally concentrated. 

sion of the analysis, taking a more consistent account of the 
surface layers of the star, augmenting the number of diag- 
nostic frequency combinations used in the calibration, and 
adding a second starting reference solar model to demon- 
strate the insensitivity of the iterated solution to starting 
conditions. We fit our model of th e frequ encies to the Bi- 
SON data discussed by iBasu et al] |20o3): they are mean 
frequencies obtained over the 4752 days from 1992 Decem- 
ber 31 to 2006 January 3 of modes of degree / = 0, 1, 2, and 
3, adjusted to take some account of solar-cycle variation 



2 THE CALIBRATION PROCEDURE 
2.1 Introductory remarks 

Naively fitting eigenfrequencies of parametrized so- 
lar models to observed solar oscillation frequencies is 
temptingly straightforward, and was one of the earli- 
est procedures to be adopted in t he present context 
IChristensen-Dalsgaard fc Goughlll98ll ). However, it is un- 
wise to adopt so crude a strategy because the raw fre- 
quencies are affected by properties of the Sun that are 
not directly pertinent to the particular inv estigation in 
hand, as w as quick ly realized at the tim e (e.g. iGoughlFigSSl : 
IChristensen -Dalsg aard fc Goughlll984h . An example is the 
effect of the near-surface layers, unwanted here, yet a se- 
rious contaminant because the region is one of low sound 
speed. It is more prudent to design seismic diagnostics that 
are sensitive only to salient properties of the structure. This 
we accomplish by noticing the roles of various structural fea- 
tures in asymptotic analysis, and relating functionals arising 
in that analysis to corresponding combinations (not neces- 
sarily linear) of oscillation frequencies. It is these combina- 
tions that are then used for the calibration. 

We emphasize that the calibration is carried out by pro- 
cessing numerically computed eigenfrequency diagnostics in 
precisely the same manner as the observed frequencies. Af- 
ter the diagnostics have been designed, asymptotics play no 
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further role. The precision of the calibration itself is inde- 
pendent of the accuracy of the asymptotic analysis; it is 
only the accuracy of the conclusions drawn from these cal- 
ibrations that is so reliant, for those conclusions depend in 
part on the degree to which the diagnostic quantities of, in 
our present study, age and heavy-element abundance, are 
divorced from extraneous influences. 



2.2 Diagnosis of the smoothed structure 

The principal age-sensitive diagnostics are contained in the 
asymptotic expression 



(71 + I L + e)vo — 
FL'^-GL'^+HL^ 
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in which i — (n, I) labels the mode, L = 1 + 1/2, and the coef- 
ficients ;= (I'o, e. A, B, I), P = 1, 11, are functionals 
of the solar structure alone, independent of i. This formula 
can be obtained by expanding in inverse powers of frequency 
the coupled pair of second-order differential equations gov- 
erning the linearized a diabatic o scillat ions of a spherically 
symmetric star, as did iTassoull (|l980l ). and at each order 
solvi ng the resulting equation-pairs successively in JWKB 
l|Gou gh 2007) approximation. Alternatively, perhaps more 
conveniently, but maybe less accurately, one can adopt an 
approximate second-order equation (which takes into ac- 
count the perturbed gravitational potential only p artially) 
and expand it alone in the limit n/L — >■ oo (e.g. iGoughl 
[l98_6b, 199^ ). The formula ((l| approximates the actual (adi- 
abatic) eigenfrequencies, for finite n, only if the scale H of 
variation of the background equilibrium state is everywhere 
much greater than the inverse vertical wavenumber of the 
oscillation mode. That is accomplished by regarding the so- 
lar model, M, to have been replaced by a smooth model. 
Ma, from which the acoustic glitches have been removed. 
We denote its frequencies by Vsi- 

The coefficients in expression ([1} that are most sensitive 
to the stratification of the core are those multiplying the 
highest powers of L at each order in uo/vsi, namely A, C, 
and F. (The L-dependent part of the leading term is also 
sensitive to the core, but merely to indicate, in the spherical 
environment, that there is no seismically detectable physical 
singularity at the centre of the star; there is, of course, a 
coordinate singularity in spherical polar coordinates.) The 
next terms in core sensitivity are D and G, and then H. 
These are also sensitive to the structure of the envelope, so 
we ignore them in the calibration. Below the near-surface 
layers of a spherically symmetrical star the integrands for 
A, C and F (which here we denote by the parameter a — 
1, 2, 3 respectively) are given approximately by 
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llGoughll201ll l. where r is a radial co-ordinate and c is the 
adiabatic sound speed; they are plotted in Fig.[T] Notice that 
the higher the order in the expansion, the more concentrated 
near the centre of the star is the integrand of the most sen- 
sitive functional. The integrands depend on progressively 
higher derivatives of the sound speed. Moreover their evalu- 
ation by fitting formula (1) to oscillation frequencies is more 



susceptible to frequency errors. Granted that we use frequen- 
cies of modes of only four different degrees, Z=0, 1, 2 and 
3, we cannot even in principle determine from them coeffi- 
cients arising in terms of higher order than those presented 
in the truncated expansion ([l]). 

One can see from expression ((2]) for the integrands of the 
coefficients A, C, and F that they depend also on vo, which 
is sensitive to the outer layers of the star where the sound 
speed is low. We remove that sensitivity by eliminating uo 
from expression ([2]) , and using instead for our diagnostics the 
parameters A — uqA, C — VqC and F = v^F, which are 
the natural factors arising in the asymptotic expansion {TJ 
of Vsi in inverse powers of Vsi/^o- 



2.3 Glitch contributions 

The abrupt variation in the stratification of a star (relative 
to the scale of the inverse radial wavenumber of a seismic 
mode of oscillation), associated with the depression in the 
first adiabatic exponent 71 = {d\np/ d\n p)s (where p, p and 
s are pressure, density and specific entropy) caused by he- 
lium ionization, imparts a glitch in the sound speed c(r), 
which induces an oscillatory component in the sp acing of 
the ei genfrequencies of low-degree seismic modes (|Goughl 
Il990al l. The amplitude of the oscillations is an increasing 
function of the helium abundance Y, and, for a given adi- 
abatic 'constant' n/p^^, is very nearly proportional to it 
(jHoudek fc Goughll2007b[ ). It is therefore a good diagnostic 
of y. To determine the amplitude we construct a deviant 



5ui 



(3) 



from the frequency Usi of a similar smoothly stratified star, 
presuming that Usi is described approximately by equa- 
tion HJ. 

A convenient and easily executed procedure for estimat- 
ing the amplitude of the oscillatory component is via the 
second multiplet-frequency difference with respect to order 
n amongst modes of like degree I: 



- 2u„,i + v.. 



n + l,i ■ 



(4) 



Taking such a difference suppresses smoothly varying (with 
respect to n) components. The oscillatory component in 
^iv, produced by an acoustic glitch, has a 'cyclic frequency' 
approximately equal to twice the acoustic depth 



dr 



(5) 



of the glitch. The amplitude depends on the amplitude F 
and radial extent A of the glitch, and decays with v once the 
inverse radial wavenumber of the mode becomes comparable 
with or less than A. 

The effects on the frequencies of a solar model M of 
a specific glitch perturbation ^71 can most readily be esti- 
mated from a variational principle in th e form v = /C/J , 



mated i rom a variational pri 
as hav e [Gqugiil (Il990al). IHou 
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Houd ck (2004 



2004|) and lMonteiro fc ThompsonI (|2005l ). 



Houdek fc Gough 



Houdek fc Gough 



(j2007b|) have found that a good approximation to the out- 
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Figure 2. The symbols in the upper panels denote seeond 
differenees A2iiy := i^n— 1 I ~ 2i'n I + ^n+l I oi low-degree modes 
obtained from the BiSON (Basu et al. 2007). The solid curve is 
a fit of the seismic diagnostic (equation ll7l l to the data by appro- 
priately weighted least squares. The dashed curve is the smooth 
contribution, including a third-order polynomial in 1/7^ to repre- 
sent the upper-glitch contribution from near-surface effects and 
a contribution from the (leading-order) second differences of 
given by equation {TJ, as descibed in the text (§2.3). The lower 
panel displays the remaining individual oscillatory contributions 
(with zero means) from the acoustic glitches to A2iU: the dotted 
and solid curves are the contributions from the first and second 
stages of helium ionization, and the dot-dashed curve is the con- 
tribution from the acoustic glitch at the base of the convective 
envelope. 



(7) 



(8) 



I -.— I p r dr 

is the mode inertia and 
6^K, ~ J S'yi p(div^)^r^ dr . 

The function ^ is the displacement eigenfunction associated 
with either or a corresponding smooth model; here we 
implicitly use Ma- Several terms in equations ([6]) and ((8]) 
are missing from the exactly perturbed equation; these are 
relatively small, and in any case to a substantial degree they 
cancel. 

The next step of the estimation is to select a conve- 
nient representation for Syi. Several formulae hav e been 
sugge s ted and used, by e.g., Monteiro fc Thompson! |l998l. 



iBasu et al. asu fcMandell (|2004l ). and lGoughl 



not all of which are derived directly from e xplicit 
acous tic glitches representing helium ionization (e.g. iBastJ 
1 19971 ). Gough used a single Gaussian function; in contrast, 
Monteiro & Thompson assumed a t riangular form; B asu et 
al. adopted a simple discontinuity (|Basu et al.lll993 ). The 
artificial discontinuities in the sound speed and its deriva- 
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Figure 3. The symbols denote contributions Su^ to the frequen- 
cies Uj produced by the acoustic glitches of the Sun (see also 
iHoudek fc GoughllioOQbh . 



tives that the latter two possess cause the amplitude of the 
oscillatory signal to decay with frequency too gradually, al- 
though that deficiency may not be immediately noticeable 
within the limited frequency range in which adequate as- 
teroseismic data are or will imminently be available. The 
analytic rep resentation, n amel y a Gaussian fu nction, which 
was used bv lGouehl (|2002l ) and lHoude^ (|2004l ). can be made 
to fit the glitch frequency perturbation more closely, espe- 
cially if the frequency range is large. 

All these early representations addressed only 
the second stag e of h elium ionizat ion. Subsequently 
iHoudek fc Goughl (|2004 |2006| . l2007al lbl) added another 
Gaussian function to take account of the first stage of 
helium ionization, relating its location, ri, amplitude factor, 
Fi, and width, Ai, to those of the second stage according 
to a standard solar model; and thereby they attained 
considerable improvement. Accordingly, we adopt that 
procedure here, and set 



71 



-{T-Tiy /2Af 



2tv ^ Ai 

1 — 1 



(9) 



summing over the two stages i (= I and II) of ionization. We 
set FiAi/riiAii = /?, n/ni — fj, and Ai/An ~ fj,- We have 
found that /3, D and fi hardly vary as Yq and Iq are varied 
in calibrated solar models, and we set their values to be 
the constant valu es 0.45, 0.70, and 0.90 res pectively, which 
gives the best fit (|Houdek fc Goughl [ioOTbl ). The quantities 
Tii, Fii and An, or equivalently ri, Fi and Ai, are adjustable 

parameters of the calibration. 

Following iHoudek fc Goughl (|2007bl ) we estimate the 
components of the displacement eigenfunction ^ of a mode 
of oscillation of Ms, and the divergence, in separated form 
as products of spherical harmonics and functions of radius 
r, using the (hybrid) JWKB asymptotic approximation (e.g. 



mg 1 

lGoughll2007h for high order n 

1/2 



cost/; , 



1/2 



Ai(-a;),(10) 



where ^(r) is the r-dependent factor in the vertical com- 
ponent of ^, having effective vertical wavenumber K, and 
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Lo = 27r!/ is the angular frequency of oscillation; the argu- 
ment X of the Airy function Ai is given by 

I 3 12/3 

X ■- sgn(?/))|-?/>| 

in terms of the phase ipir) = f K dr, which we approximate 
using a plane-parallel polytropic envelope of index m: 



(11) 



fcwf - (m-f l)cos-i (^) forf>rt, 
\K.\u}f - (m + l)\n{^ + \k.\) forf<rt, 



(12) 



in which f = t + cj^^e, with e being a phase constant, and 
Tt is the associated acoustical depth of the upper turning 
point, at which the wavenumber K vanishes. The function 



1 - 



m + 1 



1/2 



(13) 



results from approximating A' as c ^{uj^ — oj^)^^'^ in which 
the acoustic cut off frequency ujc is approx imated by (m -f 
l)/f. Following iHoudek fc Goughl (|2007bh we take m = 
3.5. The Airy function must be adopted in the expression 
for div^, which appears in the integral for S^ylC in equa- 
tion (|10p . because the upper turning point of the highest- 
frequency modes is within the He I ionization zone where 
^71 is nonzero. It is adequate to use the sinusoidal (JWKB) 
expression for both ^ and the horizontal component of the 
displacement ^ - which is determined as a horizontal deriva- 
tive in div^ - in computing the inertia, given by equa- 
tion (I7|), because almost all of the integral comes from re- 
gions far from the turni ng points. It is approxim ated by 
1 ~ ^Tui - i(m + l)n jHoudek fc Goughl boOTbl ). where 
T = r(0) is the acoustic radius of the star. The phase factor 
e was introduced to take some account of the variation with 
u} of the location of the upper turning point. 

Inserting these expressions into equations ©-(HI) yields 
the following approximation to the helium-glitch frequency 
component: 



27rylnA„^ [i^+^{m+l)uo] 

T 

^.-1 -(-^-^ni)^/2M^A'* 


T 

«:ri'e-^"-""''/'^"|2;|5jAi(-a;)i^dr 





'"Ixj^ jAi(-x-)rdr 



(14) 



where Ki :— t^iji), and where we have introduced a frequency 
amplitude factor An — iPuT"^. 

There are three additional components to l^2Vi that we 
must consider. The first is due to the abrupt variation in 
the vicinity of the base of the convection zone at Tc. We 
model it with a discontinuity in cj^ at Tc coupled with an 
exponential relaxation to the smooth model Ms vci the ra- 
dia tive zone beneath, with ac oustical scale time to= 80s, as 
did lHoudek fc Goughl (|2007bD . This leads to 



5cV 



A,ulv~^ (1 + l/lG-K^rW) 



-1/2 



X jcospi/'c +tan"^(47rro!/)]-(167r^fi'i^^-H)^''^|(15) 

where i/jc '/'(^c) and fc ~ '?(''"c), and Ac is proportional 
to the jump in uj^. 

The other two components, whose sum we denote by 
Su^i, contain a part that is generated in the very outer lay- 
ers of the star - by the ionization of hydrogen, the abrupt 



stratification of the upper superadiabatic boundary layer 
of the convection zone, and by nonadiabatic processes and 
Reynolds-stress perturbations associ ated with the oscilla- 
tions, which a re difficult to model fe.g. lRosenthal et al.lll995l : 
iHoudekllioiol l - and a part that results from the incomplete 
removal of the smooth component when taking a second dif- 
ference. The latter was obtained from equation U]), and is 
given approximately by the second derivative of v^i with 
respect to n, regarded as a continuous variable, retaining 
only the leading term. The degree-dependent term is much 
smaller than the other, and it is adequate here to regard 
the entire contribution as part of the essentially degree- 
independent upper (near-surface) glitch term, even though 
it actually arises in part from refraction in the radiative in- 
terior. We approximate it as a series in inverse powers of v, 
truncated at the cubic order: 



(16) 



We appreciate that in principle there should be an addi- 
tional contribution from the stellar atmosphere which, be- 
cause it is produced far in t he upper evanescent region of the 
mode , is a high power of v (jChristensen-Dalsgaard fc Goughl 
Il980l ). However, for the Sun and Sun-like stars its con- 
tribution to the second differences, used for determin- 
ing Vj, is small, as can be adduced from the work by 
iKield sen. Bedding fc Christensen-DalsgaardI (|2008l ). Its ef- 
fect on the fitting of the smooth components v^i is mainly 
to distort the values of B, E and I. However, these co- 
efficients are not used for the Iq and Yo[Zo) calibration. 
Accordingly, we can safely ignore this surface con tribution. 
iDogan. Bonanno fc Christensen-DalsgaardI l|201ll ) have re- 
cently illustrated this general point with specific numerical 
examples. 

The complete expression for the second difference 



(17) 



was then fitted to the second differences of the solar, or 
solar-model, frequencies to determine the coefficients rja = 
(ylii, An, m, eii, Ac, Tc, Ec, ao, ai, ^2, as), a = 1, 11. 
From the outcome, putative frequency contributions Sui'i 
were obtained by summing the second differences H16p to 
yield 

SuUi ~ A + Bvi 



X {Vn + 1,0 — Vn-l,o) 



(18) 



The initially arbitrary constants of summation A and B 
were selected in such a way as to minimi ze the L2 norm of 
SuUj, n amely ^^{A + Bui + F^i)^, as did lHoudek fc Goughl 
(|2009ah . 

The fitting of the second-differences was accomplished 
by minimizing 



(19) 



using the value of vo obtained from the fitting of expres- 
sion (|l|. (That fitting was accomplished by minimizing the 
appropriately weighted mean-square difference Es - defined 
in § 2.4 - from the smooth frequencies i/gi, which are them- 
selves derived from the raw frequencies by subtracting the 



© 0000 RAS, MNRAS 000, 000-000 



6 G. Houdek & D. 0. Gough 



A 

10-- 
9- 

I \ \ ^ 

1 2 3 4 age 

8- 

7- 

6 - 

5- 

Figure 4. Denotation of the eleven solar models which wo 
have used for calculating the partial derivatives H^p, calibrated 
to a present radius Rq = 6.9599 X 10^" cm and luminosity 
L© = 3.846 X lO^^ergs"^. The 'central model' is Model 0; the se- 
quence of the five models (0-4) has a constant value of Zo = 0.02 
but varying ago U (4.15, 4.37, 4.60, 4.84, 5.10) Gy; the second 
sequence, of seven models (0,5-10), has constant ago i* = 4.60 Gy 
but varying Zq (from 0.016 to 0.022 in uniform steps of 0.001). 
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Figure 5. Integrand — (dc/da;)/x of A as a function of radius 
fraction x := r/R for the calibrated solar models 5, and 10 
with varying Zg at constant age tQ = 4.60 Gy. The dashed, solid 
and dot-dashed curves are the results for models 5 (Zo=0.016), 
{Zo = 0.020) and 10 (Zo=0.022) respectively Note that the 
sensitivity to Zq lies mainly near the centre; the same is so for 
the sensitivity to U llGoueh fc NovotnvllTgga) . 



glitch contribution obtained by minimizing ; the two mini- 
mizations were carried out iteratively in tandem.) Here C^^^- 
is the element of the inverse of the covariance matrix 

Ca of the observational errors in A2ii', computed, perforce, 
under the assumption that the errors in the frequency data 
Ui are independent. The resulting covariance matrix Cna-y 
of the errors in ria was established by Monte Carlo simula- 
tion, using 6000 realizations of Gaussian-distributed errors 
in the raw data with variance in accord with the published 
standard errors. In carrying out the simulations we omitted 
the surface term 5u!^i, which has insignificant influence on 
the statistics. 

The outcome of the fltting to the BiSON data is dis- 
played in Fig.O the upper panel displays the second differ- 
ences, together with the complete fitted formula H17p (solid 
curve) and its individual smooth frequency contribution Sui'i 
estimated by equation (I16p : the corresponding oscillatory 
frequency contributions A2S^i' (dotted and solid curves for 
the two stages of helium ionization) and A25cV (dot-dashed 
curve) are illustrated in the lower panels of Fig.[51 All the 
frequencies displayed in the figure have been used in equa- 
tion (|19|) for fitting expression (|17|) . 

In Fig. [3] is displayed the sum of all acoustic glitch 
contributions to the frequencies estimated by fitting equa- 
tion (1 1711 to the low-de gree solar frequencies observed by 
BiSON lsasuet al.ll2007l '). 

2.4 Calibration for age and chemical composition 

We subtract the glitch contributions Sui from the full fre- 
quencies to obtain corresponding glitch-free frequencies Vsi- 
The procedure is carried out for the solar observations, for 
the eigenfrequencies of the reference solar model, and for 
the grid of models used for evaluating derivatives of the fit- 
ting parameters with respect to t* and Zo (see Fig.[4|. Then 
we iterate the parameters defining the reference model by 
minimizing Eb := (i^si — 'S'i)C^^ (I'sj — Sj), where Ca is the 



covariance matrix of the statistical errors in i/^i, which are 
determined from the independent observational errors in Ui 
and the covariance matrix C,,^^ , to obtain both the coeffi- 
cients and the covariance matrix C^0s of their errors. In 
this iteration process for only glitch-free frequencies fsi 
with k = n + > 15 were considered, because the asymp- 
totic expression ([T]) is not sufficiently accurate for lower val- 
ues of k. Each component of ^/j is an integral of a function 
of the equilibrium stratification. Some of these are displayed 
in Fig.[T] The integrals A, C and F are those of particular 
importance to our analysis, because C and F are dominated 
by conditions in the core, and, although the contributions 
to A from the core and the rest of the star are roughly 
equal in magnitude (and potentially have opposite signs), 
th e contribution from the en velope is relatively insensitive to 
U (|Gough fc Novotnvlll990l ) and Zo (Fig.p. The integrands 
in the remaining integrals are either more evenly distributed 
throughout the Sun or are concentrated near the surface. 

The differences between the smoothed frequencies v^i 
and the fitted asymptotic expression Si given by equa- 
tion lU are displayed in Fig. [S] for the BiSON data (left 
panel) and for the central model mO (right panel). 

We have carried out age calibrations using various com- 
binations of the parameters 

Cc = (A,(;,F,-57i/7i), a = l,...,4, (20) 

where —571/71 = ^ii/\/27rfoAii is a measure of the maxi- 
mum depression in 71 in the second helium ionization zone, 
and which for convenience we sometimes denote by F. The 
values of — <57i/7i and the asymptotic coefficients A,C,F 
appearing in expression ((l|, determined from the observed 
seismic frequencies, are listed in Table [T] for the Sun, and 
are plotted in Fig. [7] for the eleven calibrated grid models. 
Presuming, as is normal, that the reference model is para- 
metrically close to the Sun, we first carry out a single itera- 
tion by approximating the reference value by a two-term 
Taylor expansion about the value (^® of the Sun: 
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Figure 6. Differences between smootlied frequencies v^i, of the Sun (left) and central Model (right), and the correspondingly fitted 
asymptotic expression Si given by {ij. Modes of like degree I are connected by solid lines. Contrary to the opinions of some commentators, 
the curves are not like the oscillatory contribution from 571 /71 , for rather than decreasing, the amplitude of the undulation increases 
with increasing v (cf. Fig. 2), indicating discrepancies in the very outer layers. 



Table 1. Asymptotic coefficients for the Sun obtained from (lin- 
ear) fitting to the glitch-free (smoothed) BiSON frequencies Vs 
the expression JTJ by weighted least squares Eg. 



UQ (/iHz) 


A 


C 


F -571/71 


136.71 


0.3005 


1.912 


69.83 0.04538 


.0 ( dCc 


) At* 







where A and IS.Z0 are the deviations of the age Iq and 
initial heavy-element abundance Zq of the Sun from the cor- 
responding values of the reference model; e^a are the formal 
errors in the calibration parameters, whose covariance ma- 
trix C^aii can be derived from C^^^ and Cja^.. A (para- 
metrically local) maximum-likelihood fit then leads to the 
following set of linear equations: 

HajCQ^pHpkQok ~ HcijC^^ijAo0 , (22) 

in which Qk = (At*, A^) + eek = Qok + eefe, k = 1,2, is 
the solution vector subject to (correlated) errors eefe, and 
A/3 = — C^j + ^C/S = ^0/3 + fC/9i the partial derivatives are 
denoted by H^j = [{dCc/dU)zo, (c>Cc/c>^)t J, j = 1,2. 

A similar set of equations is obtained for the formal 
errors eek- 

HajC(^^pHpkeBk = HajC^^^peQis , (23) 

from whose solution the error covariance matrix Cekq = 
eekseq can be computed. 

The partial derivatives Haj were obtained from the 
set of eleven calibrated evolutionary models (see Fig.|4]) 
of the Sun that were u sed in a similar calibration by 
iHoudek fc Gouehl i 2007al ) . The models were computed 
with t he evolutionary programme bv lChristensen-DalsgaardI 
(|2008h . adopting the Livermore equation of state and the 



OPAL92 opacities. The set comprises two sequences: one has 
a constant value of the heavy-element abundance Zo = 0.020 
but varying age (t* = 4.15, 5.10 Gy in uniform steps of 
0.05 in Int*); the other has constant age t* — 4.60 Gy but 
varying Zo (Zo=0.016,.. .,0.022 in uniform steps of 0.001). 
Note that, for prescribed relative abundances of heavy ele- 
ments, the condition that the luminosity and radius of the 
Sun agree with observation defines a functional relation be- 
tween Yo, Zo and t* amongst the models. In Fig.[7]are plotted 
the seismic parameters A,C,F and —571/71 of the eleven 
models, each calibrated to the solar radius and luminosity, 
for determining the partial derivatives Haj of A, C, F and 
—571/71 with respect to stellar age t* and initial heavy- 
element abundance Zo- The values of the partial (logarith- 
mic) derivatives Haj so obtained are listed in Table[2] No- 
tice that within the range of model parameters that we have 
considered, the derivatives are almost constant. 



3 RESULTS 

Provided that the reference model is close to the Sun, the 
single iteration described in the previous section should pro- 
vide as reliable an estimate of (At0, AZo) as the calibra- 
tion is currently able to provide. We therefore discuss at 
first the results of single iterations. Calibrations were car- 
ried out using different combinations of the parameters C^a 
and two different reference models. They are summarized in 
TableO The older reference model is the central 'Model 0' 
which has age t* = 4.600 Gy; the second is 'Model 2', which 
has an age = 4.370 Gy. Because the acoustic properties 
of the stars in the grid vary almost linearly with t* and Zq 
the error covariance matrices associated with the single iter- 
ation are in distinguishable. We adopted the sa me physics as 
in Model S (jChristensen-Dalsgaard et al.lll99^ ) in the evolu- 
tionary calculations of both models. We notice by comparing 
rows 4 and 6 with rows 5 and 7 in Table|3]that calibrations 
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Figure 7. Seismic parameters vq, A,C and — 571/71 of eleven calibrated solar models from which their partial derivatives H^j v^ith 
respect to stellar age and initial heavy-element abundance Zq, listed in Table [21 are obtained. 

Table 2. Partial logarithmic derivatives {Haj) obtained from the two sets of calibrated evolutionary models for the Sun dcnotatcd in 
Fig. 4. They were determined by linear regression about the central Model 0. 
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without 571/71 are less stable to a change in the reference 
model than are the calibrations including 571/71. They are 
possibly less reliable, for the reasons explained in the in- 
troduction, although the result may perhaps be simply a 
symptom of slower convergence. 

To ascertain whether the entire calibration procedure 
converges, we have performed several additional iterations. 
At each iteration the corrections At, and AZq are used to 
define parameters of a new reference model, which is then 
constructed by performing another evolutionary calculation, 
followed by the evaluation of a new set of corrections At, 
and AZo as before. We repeated this for five iterations, for 
each of the two reference models, obtaining the two 'final' 
reference models, listed in Table [S] for two different com- 
binations of ("d, whose present heavy-element abundance 
Z is displayed in Fig. 1111 The progressive corrections At* 
and AZo are plotted in Fig.[Sl In carrying out the itera- 
tions we did not recompute the partial derivatives Haj and 
the corresponding error covariance matrices. To have done 
so would have been computationally much more expensive, 
would have been likely not to have speeded up convergence 
by very much, and would not have altered the final solution. 



The final residuals As :— Vsi — Si are barely distinguishable 
from those from the Sun, as illustrated in Fig.lHl 

Error contours corresponding to the calibration from 
Model in the first row of Table[S] are plotted in Fig. llOl 
Corresponding contours for Model 2 are the same, except 
that their centres are displaced to (4.603 Gy, 0.0155). One 
can adduce from our description of the analysis in ij |2.4l that 
our current treatment of the errors is not completely unbi- 
assed, because, aside from i/q, we assess the error covariances 
of the parameters defining the smooth and the glitch com- 
ponents independently; however, the potential bias is of the 
order of only {Svi/vil or less, which is small. 

Fig-El] depicts the heavy-element profiles after five it- 
erations from the two reference models (Models and 2). 
Both models have a surface value Za = 0.0142 ± 0.0005, 
which is a bout 6% higher than the value of Zs =0.0134 re- 
ported by lAsplund et al.] (|2009l ) and about 9% smaller than 
the va lue of Z^ = 0.0156 ± 0.0011 reported by ICaffau et"al] 
liooi). The error-bars of Caffau's Za value, obtained from 
numerical simulations, is indicated by the shaded region. 

The calibrated age inferred from Model after five 
iterations is 4.604±0.039 Gy, and that from Model 2 is 
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Table 3. Age calibrations with different combinations of C,^ and for tlie two reference models: Model with an age tj, = 4.60 Gy and 
Mode l 2 with an age t* = 4.37 Gy, both computed with an initial heavy-element abundance Zq = 0.02 adopting the lGrevesse fc Noelsl 
lll993l ) solar composition. The first three columns show the results adopting Model as the reference model, the fourth, fifth and sixth 
columns display the results for Model 2. The remaining three columns are the values of the error covariance matrix Ce for both reference 
models. 
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Figure 8. Corrections of the solar age tQ (left panel) and initial heavy-element abundance Zq (right panel) for six calibration iterations 
and the combination A,C,F and —57/7 for two reference models with an initial age of t = 4.60 (solid black curves) and t = 4.37 Gy 
(dashed red curves) and initial heavy-element abundance Zq = 0.02. The horizontal dashed lines are the estimated 1 cr error bars of the 
calibrated age and initial heavy-element abundance. 



4.603±0.039 Gy, using the parameter combination A, C , F 
and —571/71. The corresponding calibrations from Models 
and 2 for the combination C, F and — 571/71 are 4.602 ± 
0.044 Gy and 4.601 ± 0.044 Gy, respectively. Tabled sum- 
marizes the calibrations after five iterations from reference 
Models and 2. 



4 DISCUSSION 

In attempting to estimate the main-sequence age of the Sun 
it is prudent to adopt diagnostic quantities that are insen- 
sitive to properties that one believes not to be directly per- 
tinent. As the Sun evolves on the main sequence it converts 
hydrogen into helium in the core. According to theoretical 
models it liberates energy at a rate whose dependence on 
time, measured in units of the age t^, is not very sensitive to 
uncertain parameters defining those models, such as the ini- 



tial heavy-element abundance Zo, provided that the models 
have been calibrated to reproduce the luminosity and radius 
observed today (cf. Gough, 1990b). The same is true of the 
quantity of hydrogen consumed, mainly because the nuclear 
relations are dominated by a single branch of the pp chain, 
namely ppl, for which there is a tight link between fuel con- 
sumption and thermal energy release. Therefore one would 
expect there to be a robust link between main-sequence age 
and the total amount of hydrogen consumed: the integral 
AH ■- Jh(r)pr^dr ■- J [Xo - X (r)] pr^ dr should be 
a good indicator of the age tQ. It can be calibrated using 
seismic diagnoses of the mean molecular mass iJ,{r) provided 
that processes other than nuclear reactions that can change 
X{r), such as gravitational settling and diffusion, are taken 
adequately into account. 

In perhaps its simplest form, solar evolution involves 
computing models of constant mass in hydrostatic equilib- 
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Figure 9. Differences between residuals As := v^i — Si obtained 
from the frequencies of the Sun, Aa,0, and the calibrated model 
after five iterations from reference Model 0, As^m, where expres- 
sion Si given by Eq. JlJ was fitted to the smoothed frequencies 
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Figure 10. Error ellipses for the calibration using the combi- 
nation A,C,F, — <57i/7i and Model as the reference model: 
solutions (tQjZo) satisfying the frequency data within 1, 2 and 
3 standard errors in those data reside in the inner, intermediate 
and outer ellipses, respectively. 



rium. The models usually depend on three initial parame- 
ters: the initial abundances of, say, helium, Yo, and the heavy 
elements, Zo, and a mixing-length parameter Qc which is 
normally held constant. It is usual to fix the relative abun- 
dances of all elements other than hydrogen and helium, a 
procedure which we too have adopted here. Demanding that 
the luminosity and radius of the model agree with present- 
day observation relates two of those parameters, say Qc and 
Yq, to the third, Zq, for any t-^. Thus one obtains a two- 
parameter set of potentially acceptable models, which here 
we characterize by the values of and Zo, and which we 
attempt to calibrate with helioseismic data. 

Several diagnostics have been used in the past. As men- 
tioned in the introduction, the first to be used for a full 
calibration w ere the two mean small separations do and di 
(|Goughll200l[ ). averages over n of dn,o and d„,i, the hope 
being that the differences in the way in which the two quan- 
tities sample the core would be adequate to disentangle 
and Zq. Unfortunately, given the precision of the data at the 
time, that could not be accomplished to a useful precision. 
Moreover, by inspecting the dependence of the calibration 
on the range of frequencies over which the averages do and 
di were determined, there was evidence of contamination by 
an oscillatory component to the signatures from seismically 
abrupt variations of the stratification in and at the base of 
the convection zone. This component is particularly visible 
in second- and higher-order frequency differences with re- 
spect to n (e.g. Goughlll990al : [Ballot. Turck-Chiezefc Garcfal 
12004 iBasu fc Mandelll2004l ). 



There were several obvious improvements to the orig- 
inal calibration based solely on fitting to raw frequencies a 
smooth-asymptotic formula, such as equation ([l} or deriva- 
tives of it, that were required to be put into place in or- 
der to obtain a more reliable calibration. The first that 
we have made is to isolate much of the signal from the 
abrupt variation of the thermodynamic properties in the 
convection zone. The intention was two- fold: First, by re- 
moving the oscillatory component from the stratification one 
is left with a smooth model for which the simple asymp- 



totic expression ([T]) is more nearly valid; second, its ampli- 
tude provides an independent me asure of the helium abun - 
dance Ys in the convection zone jHoudek fc Goughll2007bl ) 
through the magnitude of the depression in 71 in the ion- 
ization zones. The latter provides, via stellar-evolution the- 
ory, the value of Yq - and therefore Xo{Zo) - in the core, 
which is required for determining the hydrogen deficit h{r). 
In carrying out the analysis, the variation in 71 has been 
represented by two Gaussia n func tions of acoustic depth, 
as rec ommended by iGoughl 1I2OO2I ) and iHoudek fc Goughl 
(|2004 . which has been found to reproduce the oscillation 
frequencies more faithfull y than eith er the simple disconti- 
nuity that was adopted by Basu ct al. (' 2004h . lBasufcMande]| 
(|200'l ) and lMazumdar fc Michel (20101'). a nd th e trian gular 
form adopted by Monteiro fc ThompsonI (|l998l . |2005| ) and 
I Verner. Chaplin fc ElsworthI ( 20061 '): presumed discontinu- 
ities in 71 or its derivatives cause the amplitude of the pre- 
dicted oscillatory feature to decay too slowly with frequency 
(IHoudek fc Goughl [2003 ). which, although apparently not 
very deleterious for the Sun, may be a serious deficiency for 
other stars. 

Another improvement is to remove from the diagnostics 
more of the infiuence of regions of the Sun that are outside 
the core. The absolute frequency of a low-degree mode 
of oscillation feels almost all of the interior structure of 
the star in inverse proportion to the sound speed along a 
ray path, except near the surface where the influence of 
the rapid variation of the acoustical cutoff frequency Uc 
dominates. The latter is largely eliminated in the small 
frequency separation, because the eigenf unctions in the very 
surface layers are almost independent of L, and therefore 
subtracting two modes of nearly the same frequency entails 
a high level of cancellation. However, the cancellation 
is not complete, simply because the frequencies of the 
two modes are not exactly the same. As luhichi (jigseT ) 
has pointed out, the ratio Ri of the small separation di 
to the large separation A; is a more direct measure of 
age, for it isolates more effective l y the nonhomologous 
aspects of the evolution (|Goughl Il990bh . and it more 
effectively eliminates the infiuence of the outermost layers 
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Figure 11. Heavy-element abundance Z as a function of tlic 
deptli-coordinate log{p) obtained from tlie reference Models 
{ti, = 4.60 Gy, solid black curve) and 2 (t* = 4.37 Gy, dashed red 
curve) after five calibration iterations. Resul ts from spectroscopic 
analyses based on nu merical simulations by I Asplund et aU ||2009| ) 
(dot-dashed line) and lCaffau et al.] l|2009l 'l (triple-dot-dashcd line) 
of the convectively driven macroscopic motion in the atmosphere 
are indicated: the s haded area indicates the reported error bars by 
ICaffau et al.l ll2009h . The initial (zero age) value Zq of both refer- 
ence models is indicated by the dotted line. After five calibration 
iterations both reference models have a surface heavy-element 
abundance Zs = 0.0142; the age obtained from the 4.60 Gy ref- 
erence model (Model 0) is 4.604±0.039 Gy, and that from the 
4.37 Gy reference model (Model 2) is 4.603±0.039 Gy (sec also 
Table Ell. 



of the Sun, as can easily be appreciated by comparing 
the formulae for d; and Ki implied by t he a symp- 
totic expression (111). i Roxburgh fc Vorontsov^ |2003l) and 
lOtf Floranes. Christenscn-Dalsgaard fc TliompsonI (|2005l ) 
have advocated that it be used for core calibration instead of 
di, ari d recently iDogan. Bonanno fc Christensen-Dalseaardl 
have illustrated its robustness numerically. Here we 
have gone further by adopting as diagnostics the factors 
A, C, and F, integrals of the solar structure which sense 
variations in conditions even more concentrated towards 
the centre of the star. 

One could consider going even further by trying to re- 
place the set of diagnostic factors with a single combina- 
tion of A, C and F designed to eliminate the influence 
of the surface layers as much as possible, analogou s to 
the p rocedure adopted bv lGough fc KosovichevI (| 19881 ) and 
iKosovichcv ct al. ( 1992.); that is tantamount to using a judi- 
ciously selected combination of small separation ratios Ri. 
Because A, C and F depend differently on the core stratifi- 
cation, the simultaneous use of all three quantities provides 
some information about the manner in which X varies with 
r. It is therefore to be hoped that the calibration is more 
secure than one using just di or Ri. It is worth mentioning 
at this juncture that the integrand for A is not actually neg- 
ligible outside the core, as can be seen from Fig. [T] indeed it 
has been known for some time that the integrand continues 
to the surface with approximately the same magnitude as it 
has at r/R = 0.5 (|Goughlll986bl see also Fig.[5]), and that 



and Zo vary, and therefore that A is at least a fairly good 
diagnostic for our purposes. We note, however, that there is 
some contamination from outside the core, as is hinted in 
Table [3] in which it is recorded that |Cei2| is smallest when 
A is not used in the calibration. 

It is also important to include the diagnostic V := 
—571/71, which measures the helium abundance in the 
convection zone, for that reflects a rather different aspect 
of the core structure and thereby enables a much more pre- 
cise determination oHq and Zq, as evinced by Tabled and 
which was already evident in an earlier phase of the investi- 
gation (|Houdek fc Goughll2007al '). Whether or not the out- 
come is more accurate depends on the reliability of the pro- 
cedure to account for gravitational settling, which relates 
in the surface to the value of Yq which controls conditions 
in the core. It should be pointed out also that V is not an 
uncontaminated measure of Ys, because it depends also on 
the e ntropy in the deep adiab atically stratified convection 
zone (|Houdek fc Goughll2007bl ). and perhaps in reality also 
on the existence of an intense magnetic field (see below). 
Our procedure could be made more reliable if we could find 
an alternative diagnostic that senses Ys more directly. 

Further remarks about the influence of the outer lay- 
ers, or the elimination thereof, are in order: In fltting the 
resolvable glitch contribution to the data an approxima- 
tion to the unresolvable contribution from hydrogen ion- 
ization and the upper superadiabatically stratified bound- 
ary layer was included, equivale nt to a cubic form in 
added to the second differences (|Houdek fc Goughll2007bh . 
Associated with the resolvable glitches are smooth contri- 
butions which were ignored in th e initial calibration for Iq 
and Zo (|Houdek fc Goughll2007al ). Subsequently they were 
taken expli citly into account, thereby removing a bias in the 
procedure (|Houdek fc Goughll2009a| [bh and, it is to be sup- 
posed, improving the accuracy of the calibration. It should 
be mentioned, however, that we have not taken explicit ac- 
count of putative errors in our modelling of the out e rmost 
layers of the Sun. IChristensen-Dalsgaard fc Goughl (|l980l ) 
found that were the oscillations to be adiabatic, the ef- 
fect of the atmosphere would be to add to the frequen- 
cies a term 5 that is itself a rapid function of frequency; 
S (yi v'' with h — 2(m + 1) for u/uc <C 1, where Uc is the 
cyclic cutoff frequency Ltj^jl-K and m is an effective poly- 
tropic index in the vicinity of the upper turning point, which, 
from fitting a (smooth) asymptotic frequency formul a to so- 
lar da t a, is expected to have a value of about 3 (| Goughl 
Il986al ; iHoudek fc Goughl l2007bl ): furthermore dln<5/dlni/ 
decrea ses with increasing v as v jvc approaches and e xceed s 
unity. iKield sen. Bedding fc Christensen-DalsgaardI ()2008h 
found that in the Sun h decreases to about 4.9 for ~ 3 mHz, 
which is not entirely inconsistent with this finding. 

We note, furthermore, that the influence of the 
perturbations on the Reynolds stress induced by the 
oscillations also has a component that increases with 
V, but more slowly tha n the effect o f a perturbation 



the integral is dominated by conditions outside the core. 
However, it appears that only the inner parts change as 



in the atmosphere (e .g. Goughl Il986bl ; iBalmforthI 1 19921 : 
[Rosenthal eFaP 1 19951 ; iHoudekl [20101 ). This result may 
also be partially responsible for the exponent 6 being 
s omewhat less than the expected value of 2(m -I- 1). 

iDogan. Bonanno fc Christensen-DalsgaardI (|201ll ) found 
that taking the correction into account obviates the neces- 
sity to use Ri instead of d; in a simple model calibration 
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for tQ in which Zq is held fixed, and yields results similar 
to those obtained from Ri with no surface term. This sug- 
gests that our neglect of the near-surface adjustment - a 
device which we adopted to maintain a workable number 
of unknown parameters in the fitting - may not be severely 
deleterious. Nonetheless, the approximation deserves further 
scrutiny. 

We have also been somewhat cavalier in our modelling 
of the acoustic glitches at the base of the convection zone. In 
particular, we have modelled them as a simple discontinu- 
ity in the second derivative of the density together with an 
exponential recovery beneath (|Houdek fc Goughll2007bl ') to 
represent standard solar models. Again, we have taken this 
approach for our convenience; after all, the sole purpose of 
modelling the glitch was to remove it. However, we are aware 
that we have not adequately taken account of the stratifica- 
tion of the tachocline, and that by so doing we risk not hav- 
ing eliminated adequately its contribution to the frequen- 
cies, and thereby may have biassed our final result. Indeed, 
it is evident that we have not been able to fit for the rapidly 
oscillating component of the second differences to the solar 
data as well as we have to the frequencies of a standard solar 
model, suggesting that t here might be room for furthe r im- 
provement of the theory. Monteiro fc ThompsonI (|2005l ) and 
IChristensen-Dalsgaard et al. (2011) have gone some way in 
making such improvements, with the intention of studying 
the stratification at the base of the convection zone itself. 
It behoves us to do so too. In this regard we observe that 
the differences between the residuals from the Sun and the 
calibrated model, displayed in Fig.O show evidence of un- 
dersampled high-frequency oscillations that are I dependent, 
hinting that the tachocline structure might be aspherical. 

It is one of our intentions to refine our core diagnos- 
tic by combining the integrals A, C and F into a single 
quantity T which measures most closely the total hydrogen 
consumption A_ff, rather than merely using the three differ- 
ent aspects of the deficiency function h{r) in parallel. By so 
doing, properties of the core that are not direct indicators 
of age should be partially eliminated, thereby increasing the 
accuracy of the calibration; furthermore, the reduction of 
the number of final calibration parameters from four to two 
would increase the precision, although that is of secondary 
concern. The construction of the diagnostic T is a tedious, 
although, we believe, relatively straightforward task which 
we have not yet completed. 

Another of our unaccomplished intentions is to report 
on varying the bounding values fci and k2 of fc — n -\- 
between which th e mode s used in the calibration are cho- 
sen to lie, as did iGoughl (|200ll . see also [Houdek fc Goughl 
I2OO8I ). This should give a better indication of the robust- 
ness of the calibration. We have carried out a partial sur- 
vey, but we are not yet satisfied with the outcome. The 
reason is that the function defined by equation (|19p . 
when evaluated with the coefficients of a corresponding 
smooth model represented by the coefficients in the expan- 
sion ([l]), has several local minima. The calibration we re- 
port here adopts the lowest of those minima. But we have 
found that as fci and fc2 are varied the relative depths of 
the minima change, and always selecting the lowest can 
lead to sudden jumping from one to another. The situa- 
tion is super ficially not unlike the earliest direct sola r model 
calibration jChristensen-Dalsgaard fc Goughl Il98ll ) , which 



also used only low-degree modes, and for which the accept- 
able minimum had eventually to be determined from other, 
rather different, seismic data. Maybe the resolution here will 
turn out to be similar. 

The standard calibration errors quoted in Tables [3]- 
[4] and illustrated in Fig. [10] are the result of propagating 
quoted observational errors in the raw frequencies. They in- 
dicate the precision of the calibration. In the absence of in- 
formation to the contrary, we have assumed that the raw- 
frequency errors are uncorrelated. It is important to realize 
that, given that some correlation is inevitable, this assump- 
tion can not only cause the precision of the calibration to 
be overestimated, but can also lea d to bias in the results 
(|Goughlll996l : iGough fc Sekiill2002l ). The calibration errors 
evidently overestimate the precision. And, of course, they 
certainly overestimate the accuracy. 

Our calibration yields =0.0142 for the cur- 

rent surface heavy-element abundance of the Sun. 
This is significantly small e r tha n that of Model S of 
IChristensen-Dalsgaard et al.l (|l996h . which has almost the 
correct sound-speed and density distribution throughout. 
Therefore our 'best' model is wrong. This conclusion is 
borne out by the fact that the residual differences, plot- 
ted in Fig.|9l are not zero. In particular, the opacity in 
the Sun, which is what Z principally determines, and 
which is fairly reliably determined by analysis of es- 
sentially all of the seism ic modes of intermediate and 
high degree (| Goughl l2004h . is not faithfully reproduced 
by our fitted model. What does that imply about the 
values we infer for Z^ and ^0? Rather than calibrat- 
ing their models for the heavy-element abundance, others 
fe.g.JPziembowski et al.lll999l: Bonanno. Schlattl fc Paternol 



|2002| : Dogan. Bonanno fc Christensen-Dalsgaardll2011 ) have 
instead simply adopted a value for Zo or Zs that was perhaps 
acceptable by other criteria, and carried out a much more 
straightforward single-parameter calibration to estimate Iq. 
The precision of such a calibration is greater than it would 
have been had Zb (or, equivalently, Zq) been included as 
a fitting parameter, but not necessarily the accuracy, even 
if the true value of Zp ha d been adopted. This matter is 
discussed bv lGoughl (|201lh . who suggests that under condi- 
tions such as these, a simple, but admittedly not reliable, 
rule of thumb is that accuracy tends to decrease as preci- 
sion increases. One cannot be sure that that is the case here 
without a much deeper understanding of the properties of 
the models against which the Sun is calibrated. 

It is of some interest to record how the outcome of 
such single-parameter calibrations depend on the values as- 
sumed for Zq. It is summarized in Table [4] for two con- 
stant heavy-element ab undances: Zo=0. 019628, the valu e 
adopted for Model S (jChristensen-Dalsgaard et al.l Il996l). 
and Zo=0. 014864, the va lue adopted for the Asplund et al 



J2009D abundances (see IChristensen-Dalsgaard fc Houdek 
[20100 . It is evident that a lower fixed value of Zq results 



in a greater solar age: an increase of 3% associated with a 
30% decrease in Zo. This is as one would expect. Reducing 
Zo requires also a reduction of Yq at fixed age, resulting in 
a less centrally condensed star and consequently a greater 
value of A. Moreover, increasing the age at fixed Zq and Yq 
reduces A. Therefore, to maintain A constant, lowering Zq 
for the calibration must be compensated by a rise in the in- 
ferred value for Iq . It should be noted that the use of a value 



© 0000 RAS, MNRAS 000, 000-000 



On the seismic age and heavy- element abundance of the Sun 13 



Table 4. Age calibrations f or fixed initial heavy-element abundance Zq = 0.019628 (the valu e adopted for Mo del S; 
IChristensen-Da lsgaard et al.| [l 99^) and Zp = 0.014864 (the value adopted for a solar model assuming the lAsplundet"aIll2009l abun- 
dances; see also Christ cnsen-Dalsgaard fc Houdelj|2010h . Results are shown for the combinations (^a = iA,C,F), (A,C) and (A), and 
for two reference models: Model with an age t* = 4.60 Gy and Model 2 with an age t^, = 4.37 Gy. The first three columns show the 
calibrated ages and associated errors for both reference models computed with constant Zq = 0.019628; the last three columns show 
the results for the age calibration of both reference models computed with constant Zq = 0.014864. For the Zq = 0.019628 age calibration 
the reference models and derivatives Ha = (dA/dt, dC/dt, dF/dt) were obtained from a grid of 5 models with varying age adopting the 
iGrevesse fc Noelj il99^ ) solar compositi on. The Zp = 0.01486 4 age calibration used the reference models and derivatives Ha from a grid 
of 5 models of varying age assuming the lAsplundetal.1 l|2009l ) solar composition. 
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4.490 


0.061 
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of Zq that is consistent with inferences from intermediate- 
and high-degree modes is a procedure which is not available 
for calibrating stars other than the Sun. 

Of course the reliability of the results of a calibration 
can be no greater than the reliability of the models that 
are used. Thus one should address the validity of the as- 
sumptions that are made, and estimate their influence on 
the inferred values of and Zq. For making the estimate 
we note that the final calibration is based on the values of 
the parameters (a'- the coefficients A, C and F of the most 
L-sensitive terms at each order in the asymptotic expres- 
sion ([T|, and the measure V of the depression in 71 due to 
Hell ionization. Were the Sun to be spherically symmetri- 
cal and nonmagnetic, the former would be indicators of the 
acoustic stratification of the core, and the latter a (model- 
dependent) measure of the surface helium abundance which 
is related, via the solar model, to the helium abundance in 
the core in a weakly Iq- and Z-dependent way. Here we ad- 
dress the potential errors in our estimates of the values of 
these two quantities. For illustrative purposes we lump the 
first three parameters together, and consider only ^1 = A, 
together with 1^4 = V. Then, from the derivatives listed in 
Table [2] one can deduce that for small errors 5 A, 5T in A 
and r the corresponding errors in and Zo axe determined 
by 



5 In Zq 



-0.91 
-3.2 



-0.58 
3.8 



5hiA 
SlnV 



(24) 



The value of A obtained by fitting the expression fl} to 
the 'smooth' frequencies can be misinterpreted by ignoring 
asphericity, which arises principally from solar activity in 
the superficial layers of the Sun. Formally, for spherically 
symmetric stars, a change in L at fixed seismic frequency u 
is associated with a change in the depth of the lower turn- 
ing point, which is what we use to gauge the modification 
of the stratification resulting from hydrogen burning. But 
there is also a modification of the L dependence by aspheric- 
ity, which we have not taken into account here. This comes 
about because of azimuthal-order-dependent filtering pro- 
duced by whole-disc measurements does not produce mean 
multiplet frequencies but is bia ssed towards the sectoral 
modes. Here we have taken the iBasu et al] (|2007l ) correc- 



tions based on a broad single-proxy correlation, rather than 
having tried to estimate the bias from the actual latitudinal 
distribution of solar activity. However, we might look at scat- 
ter and estimate the residual. To esti mate the effect on ou r 
calibration we use the observations of IChapl in et al. I (|2007h . 
who plot mean frequency differences at different epochs, av- 
eraged over n, for different values of degree I. These can be 
fitted to to estimate the corruption S^A to the coeffi- 
cient A. Averagi ng over the interval of observations of the 
iBasu et al] (|2007l ) data set that we use for our calibration 

yields S^^A ~ —5 x 10~^, implying that Iq would be over- 
estimated by 0.06 Gy and Zq overestimated by 0.009. These 
systematic changes are not negligible: the change in Iq is 
comparable with, although somewhat larger than, the typ- 
ical random errors listed in Tables |3] and [5] in calibrations 
that use F; the change in Zq is some twenty times larger. 
Asphericity arising from the centrifugal force of rotation is 
negligible for the Sun, but it can be significant in rapidly 
rotating stars. The asphericity of the solar tachocline is also 
insignificant at our present level of precision. 

Errors in the diagnostics of Yq have two obvious 
main sources. The first is the relation between ^71 and 
Ys, which depends on the equation of state, which we 
know is not accurate to a degree of precision that we 
would like. The matter h as been discussed extensively by 
iKosovichev et aL (Il99^) . Christensen-Dalsgaard fc DappenI 



I 1992h and bv lBaturinet al.l (120001 )7 and we do not pursue 
it here. The second is the relation between Ys and Yo, and 
also Zb and Zq, which depend on gravitational settling. It 
is difficult to assess the accuracy of currently used prescrip- 
tions, and it is likely that the uncertainty will remain with 
us for some time. Yet we note that it is not unlikely that 
the uncertainty exceeds our statistical errors arising from 
data errors. We note first that Yq and Ya differ by about 
0.026. This represents the amount of gravitational settling 
out of the convection zone that has taken place over the life- 
time of the Sun. If the computation of the settling rate were 
underestimated by 20%, say - an error which we do not 
regard as being unrealistically high, given the uncertainty 
in the value of the Coulomb logarithm used in truncating 
the e lectrostatic particle interactions (cf. iMichaud fc Proffitig 
Il993h - then the effective 5V could have been underes- 
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Table 5. Age Iq and initial heavy-element abundance Zq of calibrated solar models obtained after five iterations from reference models 
Model and Model 2 using the parameter combinations = (A, C, F, —1571/71) and (C, -F, —1571/71). The values for the initial 
helium abundance Yq and surface abundances Zs and Ys are obtained from the models. The last three columns are the standard errors 
(components of the error covariance matrix C©) associated with the calibrated values of Iq and Zq. 



tQ (Gy) Zo Yo Zs -(-C©i2)i/2 c]^^^ 

Model 

A, C, F, -<57i /71 4.604 0.0155 0.250 0.0142 0.224 0.039 0.0013 0.0005 

C,F, -571/71 4.602 0.0155 0.251 0.0142 0.224 0.044 0.0004 0.0005 

Model 2 

A, C, F, -<57i /71 4.603 0.0155 0.250 0.0142 0.224 0.039 0.0013 0.0005 

C,F, -571/71 4.601 0.0155 0.251 0.0142 0.224 0.044 0.0004 0.0005 



timated likewise, and according to equation (|24|l errors of 
+0.05 Gy and -0.001 would be imparted to and Zq re- 
spectively. What is perhaps more serious is the possibihty 
of material redistribution in the energy-generating core, ei- 
ther by large-scale convection or by small-scale turbulence 
induce d possibly by rotational shear. iGough fc KosovichevI 
argued that there is evidence for that having oc- 
curred, concomitant with a reduction of the sound-speed 
gradient in the innermost regions, and thereby making the 
Sun appear younger than it really is. This matter should 
perhaps be investigated further in the future. But what 
requires serious consideration now is the degree to which 
a magnetic field might suppre ss the acoustic glitch asso- 
ciated with helium ionization. Basu fc Mandell (|2004l ) and 
IVerner. Chaplin fc ElsworthI (j2006f l have reported a diminu- 
tion during solar cycles 22 and 23 in the amplitude of 
the acoustic signature of the glitch with increasing mag- 
netic activity (gauged by the 10.7cm radio flux -F10.7), 
with an average slope dlnF/dFio.T — —0.001 (in units of 
10'^^J~^sm^ Hz). It has already been pointed out that that 
requires magnetic field strength vari ations of ord er 10^ G 
in the second helium ionization zone (|Goughll2006l '). More- 
over, it is much greater than that implied by Libbrecht & 
Woodard's (1990) observations in the previous cycle. Given 
that the average of Jio.7 over the interval of observation of 
the BiSON data was about 120, this would imply, had the 
magnetic perturbations been small, that T has been under- 
estimated by about 10%, namely 7xl0~^, implying that Iq 
has been overestimated by about 10% and, formally, Zq un- 
derestimated by about 90%. This result appears to render 
hopeless any attempt to calibra te the glitch to determine Ys. 
Howev er, we note that recently IChristensen-Dalsgaard et al.l 
l|201ll ) have found no evidence for such variation. It be- 
hoves us, therefore, urgently to investigate the matter fur- 
ther. Magnetic-field issues aside, the relation between T and 
Ya is reliant o n the equation of state, which we kiiow to be 
deficient (e.g. iKosovichevet al.|[l992l : iBaturin et al.ll2000l ). 

There are other assumptions that are implicit in most 
solar evolution calculations. Two which have obvious seri- 
ous implications regarding the apparent age are the con- 
stancy of the total mass M of the Sun - the assumption 
is that there has been no significant accretion nor mass 
loss on the main sequence - and that physics has not 
evolved such that, in appropriate units, Newton's gravita- 
tional constant G varies with time. Failure of either of those 



two assumptions can lead to a substantial deviation of the 
Sun's evolution from the usual standard. Numerical com- 
putatio ns of the effect of varying G we re carried out long 




Pochodafc SchwarzschildT 1 19641 ). lEzer fc CameronI 



Roeder fc Demarguel ( 19641 ) and IShaviv fc BahcaU 



an d computations w it h mass loss have be e n per - 
formed byTGuzik et all (Il987l). ISwensen fc Faulk ner' (l99^, 
Guzik fc Coxl (I1995I). fsackmann fc Boothroydl (, 2003) and 



Guzik fc Mussackl ()2010l )~ the resul ts are summariz ed by an 



analytical approximation given bv lGoughl lll990bl ). In par- 
ticular, if G or M had been greater in the past, then the 
solar luminosity would have been greater, and more hydro- 
gen would have been consumed, and the Sun would now be 
younger than it appears. It behoves us also to mention that 
the luminosity of the Sun today is somewhat greater than 
the standard value that we have adopted in this work, be- 
cause the surface radiant fiux increases slightly with increas- 
ing latitudinal direction. Such issues go beyond the scope of 
this investigation. 



5 CONCLUSION 

We have attempted a seismic calibration of standard solar 
models with a view to improving earlier estimates of the 
main-sequence age tQ and the initial heavy-element abun- 
dance Zo. Our long-term goal has been to achieve a preci- 
sion which could distinguish between planet formation oc- 
curring simultaneously with or subsequent to the formation 
of the Sun. Our current best estimates, around 4.60±0.04Gy, 
are summarized in Table[5] the age is close to the previ- 
ous preferred values - in particular, the age adopted for 
Christensen-Dalsgaard's Model S - and the implied present- 
day surface heavy-element abunda nce lies between the mod- 
ern spectros c opic v alues quoted bv lAsplund et al.l (|2009l ) and 
ICaffau et al.l (|2009l ). However, we emphasize that there re- 
main many uncertainties in our procedure, and that future 
revision is not unlikely. 
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